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Hydrogen peroxide (H202) has its significance during the combustion of heavy hydrocarbons in the internal 


combustion (IC) engines. Owing to its importance the measurements of H,O, dissociation rate have been reported 


mostly using the shock tube apparatus. These types of experimental measurements are although quite reliable but 


require high cost. On the other hand, numerical simulations provide low cost and reliable solutions especially us- 


ing computation fluid dynamics (CFD) software. In the current study an experimental shock tube flow is modeled 


using open access platform OpenFOAM to investigate the thermal decomposition of H203. Using two different 


convective schemes, limitedLinear and upwind, the propagation of shock wave and resultant dissociation reaction 


are simulated. The results of the simulations are compared with the experimental data. It is observed that the rate 


constant measured using the simulation data deviates from the experimental results in the low temperature range 


and approaches the experimental values as the temperature is raised. 


Keywords: Reaction kinetics, Rate constant, Numerical simulation, Shock wave, Reflected shock wave, Contact 


discontinuity, Expansion fan, Internal flow, Compressible flow, Gas dynamics 


Introduction 


During the combustion of fossil fuels in the internal 
combustion engines, numerous intermediate reactions 
occur before the production of final products (Bhaskaran 
et al. 2002). Hydrogen peroxide is one of the most im- 
portant intermediate compounds formed which guides 
and controls the downstream chain of reactions in the 
temperature range of 850 K to 1200 K (Hong et al. 2011). 
Hydrogen peroxide is also used as monopropellant using 
the catalytic decomposition technique in some types of 
rocket engines. It has its future as a green monopropel- 
lant (Westbrook 2000). For the delivery of micro satel- 
lites in the orbits, it is an established monopropellant fuel 
(Davenas et al. 2004). It is also being utilized as an oxi- 
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dizer in some rocket engines. It has also found its utility 
for the alignment of satellites and attitude control pur- 
poses. Besides hydrogen peroxide has numerous indus- 
trial applications; for instance treatment of waste water, 
industrial waste treatment, bleaching of textile products 
and bleaching of paper etc. 

Due to its versatile applications, HO, dissociate reac- 
tion is extensively used experimentally. Most experi- 
mental setups used shock tube. In the shock tube envi- 
ronment, the propagating shock wave upon reflection 
from end wall provides appropriate environment (in 
terms of temperature and pressure) where the thermal 
decomposition can occur (Anderson 2003). This process 
has been studied using various diagnostic techniques. 
Bilwakesh et al. (1968) performed the thermal decompo- 
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Nomenclature 
A speed of sound (m/s) 


C sp. heat at const. pressure (joules/kg-K) 


P 
diffusion coefficient of lth species 


specific internal energy (joules/kg) 
specific enthalpy (joules/kg) 
Identity Matrix 

diffusion flux of Ith species 
rate constant(em*mol’'s”') 
Kelvin (temperature) 
Kinetic Energy (joules) 
length (m) 

meter 

Pressure in Pascal 

Pascal (pressure) 

heat flux (Watt/m7) 


sition analysis with the help of experimental shock tube 
facility. In order to monitor the thermal decomposition 
they utilized the absorption methods. The driver section 
consisted of air and helium mixture whereas the driven 
section consisted of H2O, and nitrogen mixture. The ni- 
trogen served as the bath gas. Meyer et al. (1969) ana- 
lyzed the decomposition of N2H; and H202. They meas- 
ured the decomposition rate of HO, between tempera- 
tures 950 K to 1450 K and the pressure was kept up to 20 
bar. UV absorption was used to monitor the decomposi- 
tion rates. 

Few other studies such as Trainor et al. (1974) con- 
ducted the study of reunification of OH ion using the 
photolysis technique at low pressure and temperature and 
utilized the absorption spectroscopy. Using flash photo- 
lysis Zellner et al. (1988) conducted experiments at low 
temperature and pressure range. Forster et al. (1995) with 
the help of laser induced fluorescence monitored the ions 
at room temperature and pressure up to 150 bar. Fulle et 
al. (1996) also used the laser induced fluorescence for the 
reunification of OH radicals. The maximum pressure and 
temperature were 150 bar and 700 K respectively. 
Sangwan et al. (2012) performed the experimental study 
of hydroxyl to hydroxyl reaction with the help of UV 
absorption monitoring technique at temperatures ranging 
from 296 K to 834 K and pressures ranging from 1 bar to 
100 bar. Bahrini et al. (2012) performed the quantitative 
analysis of hydrogen peroxide. Brouwer et al. (1987) 
performed the theoretical calculations of reaction rates by 
using statistical adiabatic models. The temperature was in 
the range 200 K to 1500 K. Troe et al. (2008) used the 
ab-initio technique for the calculation of decomposition/ 
recombination of H2O, up to maximum temperature of 
5000 K. This ab-initio technique is a quantum chemistry 
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method and is a subject of computational chemistry 
(Levine 1991). Sellevag et al. (2009) used two transition 
state model for the calculation of rate constants with 
temperature ranging between 200 K and 3000 K. 

Using shock tubes, Hong et al. (2009) conducted ex- 
periments regarding the H,O, thermal decomposition. 
They used tunable laser absorption near 2.5 x 10° m to 
detect the H,O in the products. Hong et al. (2010) used 
IR absorption for the detection of H20 at 2.55 x 10° m in 
the temperature between 1000 K to 1200 K and pressure 
between 0.9 atm. to 3.2 atm. In another study, Hong et al. 
(2011) used UV absorption to monitor OH near 306.7 x 
10° m and used infrared absorption to monitor water near 
2.55 x 10°m. They measured the rate constant for hy- 
drogen peroxide thermal decomposition at temperatures 
between 1020 to 1460 K and pressure at 1.8 atm. More 
recently, Sajid et al. (2013) performed the experiments 
for the analysis of thermal decomposition of H,O, in a 
shock tube and used quantum cascade laser absorption 
near 7.7 x 10° m. They performed the experiments in 
temperature range 930 — 1235 K and for pressures at 1, 2 
and 10 atm. 

As briefed above, most of the work in this area is ex- 
perimental. On the other hand, Computational Fluid Dy- 
namics (CFD) is a fast growing area to simulate the ex- 
pensive physical experiments. Due to cheap availability 
of computers and state of the art softwares/tools which 
are both open source and freely available, it has become 
feasible to perform simulations in a rapid and cost effec- 
tive manner. To harness the potential of CFD, the expe- 
rimental setup described by Sajid et al. (2013) has been 
simulated in an open source OpenFOAM software. One 
of its solvers, reactingFoam has been used to perform the 
simulations using the limitedLinear and upwind convec- 
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tive schemes. The reactingFoam solver supports the reac- 
tion kinetics. The thermal decomposition analysis of hy- 
drogen peroxide has been conducted by modeling the 
experimental shock tube and the results have been com- 
pared with the experimental results. 


Shock Tube analytical model 


Theory 


The shock tube is a circular or rectangular pipe with 
both ends closed. The shock tube for the current study 
has circular cross-section. A thin membrane called di- 
aphragm divides the shock tube into two compartments. 
In Fig. 1 the left compartment is called the driver section 
and is denoted by ‘4’. The right compartment is called 
the driven section and is denoted by ‘1’. Both compart- 
ments are filled with gas whereas the pressure of driver 
section is greater than the driven section. Both the sec- 
tions can have the same gas or different gases. Similarly 
both compartments may have the same temperature or 
different temperature. 


Diaphragm 


4 1 


Driven Section 
P,, T), Gas 


Driver Section 
P, ‚T; ‚Gas 


P> P; 
Fig. 1 Initial Profile of Shock Tube 


The velocities in both compartments are initially zero. 
Suppose that Pı, Tı are pressure and temperature in 
compartment 1 and P,, T; are pressure and temperature in 
compartment 4. Then the pressure ratio P4/P; is said to be 
the diaphragm pressure ratio. 

Figure 2 shows different regions formed due to sudden 
removal of diaphragm by some mechanism. The region 1 
and 2 is separated by normal shock wave which is head- 
ing towards right. The region 2 and 3 are separated by 
contact discontinuity. There is an expansion fan in be- 
tween the regions 3 and 4. The regions 2 and 3 have same 
velocity and pressure but different temperature and density. 
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Fig. 2 Pre Shock Reflection Scenario 


Figure 3 is depicting the post shock reflection scene. 
When the shock is reflected back after striking the end 
wall of driven compartment, the region 5 appears. It has 
comparatively high pressure and temperature and zero 
velocity. 
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Fig. 3 Post Shock Reflection Scenario 


Governing Equations 

The movement of shock wave in the shock tube can be 
modeled analytically using continuity, momentum and 
energy balance across the normal shock wave. 


pW = p (W-v) (1) 
Pit pW =pr+ p (Wy (2) 
2 W- 2 
e ++ satel v) (3) 
P 2 P2 2 


Here p denotes density, W represents shock velocity, v 
is the velocity of contact discontinuity, p represents pre- 
ssure and e Here p denotes density, W represents shock ve- 
locity, v is the velocity of contact discontinuity, p represents 
pressure and e is the internal energy. The Eq. (4) is the Hugo- 
niot equation, whereas w represents the specific volume. 


e-e, = EP) om -w) (4) 

The Eq. (5) provides the shock strength a from 
which the unknown value of p» is obtained. | 

nfi C-D- a 


P| 2x (27 +(7+1))(p2/p1-1) Pi 


Where y is the specific heat ratio and a represents the 
speed of sound. 

The Eqs. (6) — (8) represent the mass, momentum and 
energy equations respectively for the reflected shock. 


PW, + v= pW, (6) 
Prt W, +v? =pst psW, (7) 
2 
W,+v w? 
h, U ) = hy -> (8) 


Where W,.is the velocity of the reflected shock and A 
is the enthalpy. 


Numerical modeling 


The Navier Stokes equations along with the species 
transport equation are given by 


 .¥7.(pu)= 
L +V-(pU)=0 (9) 
00) (pUU)=V-p+V-t (10) 
O( pK 
(21) .9.(pht) + E av (p ee (11) 
=-V q+V-(t U) 
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Y, 
ALP), 9. (pru)+¥-5,=R, (12) 
J, =-pDNY, (13) 
r=9[Vu+(vuy' ]-Ż89(v-v)1 (14) 


p = pRT (15) 

Where U is velocity vector, t represents viscous tensor, 

I is the identity matrix, Ois dynamic viscosity, K 

represents kinetic energy, Y; is mass fraction of /" species, 

Jıis the diffusion flux of l" species, R; represents produc- 

tion rate of the /™ species, D; is diffusion coefficient, R is 

the gas constant, q is heat flux, pis density, p is the 

pressure, and A is the enthalpy. The chemical reaction and 

the rate constant for the reaction are given by Eq. (19) 
and Eq. (20) respectively. 


Solver 


The OpenFOAM supports a diverse range of solvers 
including the pressure based and density based solvers. 
The shock tube phenomena can be modeled by sonic- 
Foam pressure based solver and rhoCentralFoam density 
based solver. But these solvers work in inert scenario 
only. As has been previously mentioned that reacting- 
Foam solver in OpenFOAM supports reaction kinetics. 
However it is a pressure based generic solver in which 
intermediate or global reactions can be modeled. There- 
fore, it has been chosen as a suitable candidate to simu- 
late the thermal decomposition of hydrogen peroxide by 
modeling the shock tube apparatus. 

The reactingFoam solver consists of PIMPLE algo- 
rithm. It is a combination of SIMPLE and PISO algo- 
rithms. The k-epsilon turbulence model has been used for 
the current study due to its robustness and accuracy. It is 
also suitable for confined flows (Versteeg 2007). In order 
to compute the dynamic viscosity, the Sutherland trans- 
port model has been used. C, is calculated using the 
JANAF tables for the species (McBride 1993). The Par- 
tially Stirred Reactor (PaSR) combustion model has been 
used to model the combustion phenomena (Correa 1993). 
It is a modification of Eddy Dissipation Concept (EDC). 
The species considered in the simulations are H202, H20, 
OH, Ar and He. 

For the time discretization, the Euler implicit scheme 
has been utilized which is first order and bounded. In 
order to interpolate the diffusion coefficient, the linear 
interpolation scheme has been utilized. To discretize the 
surface normal gradient the uncorrected scheme has been 
implemented which is first order and bounded. 


Convective Schemes 

For the aforementioned case two different convective 
schemes, LimitedLinear and upwind, are used. Limited- 
Linear falls into the helm of Total Variation Diminishing 
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(TVD) schemes (Direct 2015). It is bounded and is of 
second order. 

Figure 4 represents the control volume around P. E 
and W are east and west nodes. e is the cell interface. The 
fluxes are calculated at cell interfaces. 


Fig. 4 Control Volume around P 


Using Fig. 4, the central difference scheme is given by 
1 
be = bp + (7) (bs — br) (16) 


Where ¢,is flux at the cell interface e, g(r) is flux 
limiter and for the case of central difference g(r)=1. 
The limitedLinear scheme is obtained by replacing the 
flux limiter of central difference scheme with the Sweby 
flux limiter. The Sweby limiter is given by Eq. (17) 
(Versteeg 2007) 

g(r) =max| 0, min(1.5r,1), min (r, 1.5) | (17) 


The upwind scheme is first order bounded. It is given 
by Eq. (18) 
de = Op (18) 


Behavior of Convective Schemes 


In order to compare the results of convective schemes 
with the analytical results, each simulation for limited- 
Linear and upwind schemes has been run by using the 
reactingFoam solver. Initially laminar and inert settings 
have been taken for the simulations. Air has been taken 
as the working inert gas. The simulation is 1-D and the 
mesh resolution is 200 computational cells. 


Table 1 Shock Tube Specs 


Quantity Symbol Value 
Driver Length Lı 0.1m 
Driven Length L: 0.1m 
Total Length Li + Le 0.2m 
Driver Pressure Py 1000,000 pa. 
Driver Temperature T4 800 K 
Driven Pressure P; 100,000 pa. 
Driven Temperature Tı 300 K 


Figures 5 and 6 show almost the same results in 
both cases of limitedLinear and upwind schemes. There 
is only slight difference in the overshoot at the normal 
shock. The overshoot is more prominent in case of limi- 
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tedLinear. The normal shock wave in analytical case is 


lagging 


Both schemes have shown that the pressure and veloc- 
ity are the same in regions 2 and 3 which conform to the 
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physics of the problem. There is an over prediction in 
pressure and an under prediction in velocity in both cases. 
The contact discontinuity has been captured in both cases, 
which has been shown in the temperature and density 


Temperature (K) 


Density (Kg/m*) 
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N 
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plots. The temperature is over predicted in both regions 
but the difference is less in region 3 as compared to re- 
gion 2. The results are very good in case of density in re- 
gion 3 and there is under prediction in the case of region 2. 


Comparison of reactingfoam results with expe- 
rimental data 


This section provides the details of the results of reac- 
tingFoam simulations regarding the thermal decomposi- 
tion of hydrogen peroxide in the shock tube model. The 
reactingFoam solver has been run for both limitedLinear 
and upwind convective schemes. 


Summary of Experimental Data 


Sajid Es-Sebbar et al. (2013) experimental data has 
been taken to run the reactingFoam simulations. Sajid 
used analytical model to calculate temperature and pres- 
sure in region 5. Equation (19) represents the hydrogen 
peroxide decomposition reaction while Eq. (20) is its rate 
constant measured from the experiment. It has been en- 
tered as an input to the reactingFoam solver. 

H20; + Ar—20H+Ar (19) 
kK(T 
(7) (20) 


= 19 (16294012) 


at 1,2 atms 


xexp(—21993+301/T) (cm*mol"'s™! ) 


Simulation Results (limitedLinear) 


1-D shock tube simulations have been run using 4000, 
6000 and 8000 spatial resolutions using limitedLinear as 
the convective scheme. The purpose to use different mesh 
resolutions is to ensure the mesh resolution independence 
for thermal decom-position of hydrogen peroxide. Each 
simulation has been run for time duration of 16 millise- 
conds with the initial time step of 5 x 107 seconds. The 
shock tube specifications used are detailed in Table 2. 
The driver section comprised of helium gas. The driven 
section consists of gas mixture as mentioned in Table 3. 


Table 2 Shock Tube Specs 


Quantity Symbol Value 
Driver Length Lı 4.5m 
Driven Length L, 9m 
Total Length L, +L, 13.5m 
Driver Pressure Py 107998 Pa 
Driver Temperature T4 800 K 
Driven Pressure P; 12000 Pa 
Driven Temperature Tı 296 K 


Table 3 Mole Fraction of Gases (Driven Section) 


Gas Mole Fractions 
Argon 0.991 
H202 0.005 

HO 0.004 
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A probe has been placed at 13.495 m from the origin, 
so that the thermal decomposition data, temperature, 
pressure and velocity could be recorded. The probe loca- 
tion has been finalized to ensure the maximum tempera- 
ture of post reflection region, since the temperature falls 
drastically as the reflected normal shock gets away from 
the driven section end wall. This phenomenon of reflect- 
ed shock wave can be viewed in Fig. 7. 


H-0, Decomposition Simulation in OpenFOAM 
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The thermal decomposition profile for 250 us from the 
beginning of decomposition in the post reflection region 
5 at the specified probe for the mesh resolution 6000 is 
shown in Fig. 8. 

The velocity, temperature and pressure scenarios for 
the 250 us thermal decomposition duration have been 
shown by Figs. 9, 10 and 11 respectively. 

It has been observed that the temperature plot could 
not capture the contact discontinuity. The region 2 tem- 
perature is nearly 690 K, which agrees well with the 
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Fig. 8 H,O, mole fraction time history 
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Fig.9 Velocity profile during thermal decomposition of H202 
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temperature 696.94 K predicted by WGD calculator for 
this region. The pressure for the region 2 is 58000 Pa, 
which is comparable to 69643 Pa predicted by WGD 
calculator for this region. The pressure and temperature 
values predicted by the reactingFoam solver, with limi- 
tedLinear convective scheme, in region 5 were lower as 
compared to WGD calculator values and the experimen- 
tal values (see Table 4). The pressure and temperature 
values in region 5 as reported by Sajid et al. are also 
based on the analytical equations. The comparison of 
hydrogen peroxide mole fraction profile with the expe- 
rimental data has been analyzed in section 4.4. It also has 
been observed that the thermal decomposition profiles 
are not changing significantly by changing mesh resolu- 
tion from 4000 to 8000. So it is safe to consider the mesh 
independence of hydrogen peroxide decomposition pro- 
file at 6000 mesh resolution. 
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Fig. 11 Pressure profile during thermal decomposition of HO, 


Table 4 Temperatures and Pressures (Region 5) 


Pressure P5 


Source Temperature T; (K) 


(atmospheres) 

(Sajid et al. 2013) 1167 2.35 
(WGD 2008) 1211.45 2.3687 
limitedLinear 

4000 mesh 875.4 1.45 
limitedLinear 

6000 mesh 878.59 1.42 
limitedLinear 

8000 mesh 874.62 1.415 


Simulation Results (upwind) 


Table 5 shows the temperature and pressure values in 
region 5 produced by reactingFoam upwind case at mesh 
resolution 6000. 

It is observed that the velocity remains the same in re- 
gion 2 as has been predicted by limitedLinear based si- 
mulation. It has been observed that like the limitedLinear 
scheme, the upwind scheme also could not capture the 
contact discontinuity. The temperature in region 2 is al- 
most 700 K. The pressure in region 2 is 60000 Pa. Figure 
14 shows that the temperature for upwind case has in- 
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creased as compared to limitedLinear. It has resulted in 
more dropping of the mole fractions of hydrogen perox- 
ide in case of upwind as shown in Fig. 12. 


Table 5 Temperatures and Pressures (Region 5) 


Pressure Ps 


Source Temperature T; (K) (atmospheres) 
(Sajid et al. 2013) 1167 2.35 
(WGD 2008) 1211.45 2.3687 
upwind 
6000 mesh 911.61 1.496 
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Fig. 12 H,O, mole fraction time history (at probe location) 
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Fig. 13 Velocity profile during thermal decomposition of H203 
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Fig. 15 Pressure profile during thermal decomposition of H203 


Comparison of Results 


Figure 16 is the comparison of the hydrogen peroxide 
thermal decomposition profiles of the limitedLinear and 
upwind schemes with the Sajid et al. experimental re- 
sults. 
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Fig. 16 H,O, Decomposition Profile (at probe location) 


The Sajid et al. data is shown by ‘+’ sign, which 
serves as the reference. The results shown by limitedLi- 
near and upwind schemes are for the 6000 grid resolution. 
The plots for limitedLinear and upwind cases are for av- 
erage temperatures 878.59 K and 911.61 K respectively. 
The average pressure values in region 5 for these two 
schemes are 1.42 atm. and 1.469 atm. respectively. Sajid 
et al. reported the temperature 1167 K and pressure 2.35 
atm. for this region, which are based on analytical com- 
putations. 

Sajid et al. plot represents the exponential decay in the 
mole fraction of hydrogen peroxide, which is a 
non-linear phenomenon. Figure 16 shows that the limi- 
tedLinear and upwind cases also follow the same type of 
non-linear behavior. The limitedLinear curve overlaps the 
Sajid et al. curve from the start of decomposition point 
up to 100 s. Afterwards it deviates due to reduction in 
decay rate and finally it terminates the graph at 
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250 s to the value of 0.00168, which is higher than the 
experimental mole fraction value 0.00078 at this point. 

The upwind curve shows fast exponential decay in the 
mole fraction as compared to experimental data during 
the first half section. The second half portion of the curve 
shows the comparatively low decomposition rate. The 
curve finally ends the graph with slightly lower value 
than the experimental value of mole fraction. The fast 
exponential decay of mole fraction is attributed to the 
high temperature obtained in the upwind scheme as 
compared to the limitedLinear scheme. The latter parts of 
the curves in both cases show decrease in the decomposi- 
tion rates, which is due to sharp decrease in temperature. 
It is again mentioned here that Sajid et al. used the values 
of temperature and pressure on the basis of analytical 
equations. 

Table 6 shows the relation between rate constant k and 
1000/T, where T represents temperature. The results of 
Sajid et al. and Hong et al. are for pressure 2 atmospheres 
and 1.7 atmospheres respectively, which have been 
shown for comparison. The limitedLinear results have 
also been shown. The limitedLinear results deviate from 
the Sajid et al. and Hong et al. results at lower tempera- 
tures and approach the experimental values as the tem- 
perature is increased. 


Table 6 Rate constant vs Temperature 


Source 1000/T k 
UK) (cm`mol -s~ ) 
1.1382 1.44 x 10 
1.0763 1.95 x 108 
limitedLinear 

0.9730 4.42 x 108 
0.8857 2.95 x 108 
1.0639 1.44 x 108 
1.0029 5.16 x 10° 

Sajid et al. (2013) sai 1.61 x 10" 
0.9057 4.54 x 107 
0.8568 1.44 x 108 
0.8319 2.47 x 108 
0.9924 8.50 x 10° 
0.9812 9.77 x 10° 
0.9724 1.15 x 107 
0.9563 1.61 x 10’ 
0.9378 2.40 x 107 

Hong et al. (2011) 0.8953 5.65 x 107 
0.8913 6.77107 
0.8648 1.01 x 108 
0.8552 1.28 x 108 
0.8407 1.95 x 108 
0.8319 2.11 x 10? 
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Conclusion 


The importance of hydrogen peroxide can be gauged 
from the fact that it is produced in the bulk quantities as 
an intermediate compound during the combustion of fos- 
sil fuels in various types of engines. It controls the chain 
of intermediate reactions leading to the final formation of 
products. It is used as fuel in some kinds of rockets, 
whereas it is also used as oxidizer in some other types of 
rocket systems. It is also used in industry. 

In the present study, the shock tube model of an expe- 
rimental setup has been simulated using the reacting- 
Foam solver in the OpenFOAM. The results of limited- 
Linear and upwind convective schemes have been dis- 
cussed and compared with the experimental results. The 
thermal decomposition curves produced by the limited- 
Linear and upwind convective schemes agree well with 
the experimental curve as shown in Fig. 16. The rate 
constant k in case of limitedLinear has been computed at 
different temperatures. Table 6 shows that it deviates 
from the experimental results at low temperatures and 
agrees well to the experimental values as the temperature 
is raised. 
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